#1.0.7.10
#### import arcpy
import arcpy.mapping as mp
import pandas as pd
import sys
import shutil
import os
import numpy 
import math
reload(sys)
sys.setdefaultencoding('utf-8')

#以下是变量，每次都得改
csv_path = r"E:\yanjiusheng\word\汇报\85+18\85+18\85+18_山杨林+油松林\温度带\85+18_山杨林+油松林_3_中温带\类型csv" #类型csv路径
poly_path = r'E:\yanjiusheng\word\汇报\85+18\85+18\85+18_山杨林+油松林\温度带\85+18_山杨林+油松林_3_中温带'
#选改不是必改
target_type = []#指定类型出点、图斑，其他类型不会出
#以下是常量，只改第一次
point_temp = r"E:\yanjiusheng\xintu\chenxu\succession\total_split\temp_point.shp" #temp_point.shp 的路径
split_csv_path = r"E:\yanjiusheng\xintu\chenxu\succession\3.31" #拆分的1927个csv路径
total_split_path = r"E:\yanjiusheng\xintu\chenxu\succession\total_split\total_split.shp" #total_split.shp 的路径
total_split_temp_path = r"E:\yanjiusheng\xintu\chenxu\succession\total_split\total_temp_split.shp" #total_temp_split.shp 的路径
location_document = r"E:\yanjiusheng\xintu\chenxu\succession\loc.txt"
save_csv_path = r"E:\yanjiusheng\xintu\chenxu\succession\csv" #演替图csv路径

def got_rect_locations():
    mxd = mp.MapDocument('CURRENT')
    lsd = mp.ListDataFrames(mxd)
    data_in_layout = True
    for df in lsd :
        if not data_in_layout:
            break
        for layer in mp.ListLayers(df):
            if layer.name == "Converted_Graphics":
                data_in_layout = False
                Converted_Graphics = layer
                break

    if not data_in_layout:
        with arcpy.da.SearchCursor(Converted_Graphics , ["OID@", "SHAPE@"]) as cursor:
            rect_locations = []
            for row in cursor:
                point_list = list(set([(point.X,point.Y) for point in row[1][0]]))
                left,bottom,right,top = min([i[0] for i in point_list]), min([i[1] for i in point_list]),max([i[0] for i in point_list]), max([i[1] for i in point_list])
                rect_locations.append ((left,bottom,right,top))
    else:
        relative_locations = []
        rect_locations = []
        mimor_df = min(arcpy.mapping.ListLayoutElements(mxd, "DATAFRAME_ELEMENT"),key = lambda x:x.elementWidth)
        left,bottom,right,top = mimor_df.elementPositionX,mimor_df.elementPositionY,mimor_df.elementPositionX + mimor_df.elementWidth,mimor_df.elementPositionY + mimor_df.elementHeight

        for elm in arcpy.mapping.ListLayoutElements(mxd, "GRAPHIC_ELEMENT"): 
            if right-elm.elementWidth > elm.elementPositionX > left and top -  elm.elementHeight > elm.elementPositionY >bottom :
                relative_locations.append((elm.elementPositionX,elm.elementPositionY,elm.elementPositionX+elm.elementWidth,elm.elementPositionY+elm.elementHeight))
        minX,minY,maxX,maxY =  mimor_df.extent.XMin,mimor_df.extent.YMin,mimor_df.extent.XMax,mimor_df.extent.YMax
        for i in relative_locations:
            absolute_minX = ((i[0] - left)/(right-left))*(maxX-minX)+minX
            absolute_maxX = ((i[2] - left)/(right-left))*(maxX-minX)+minX
            absolute_minY = ((i[1] - bottom)/(top-bottom))*(maxY-minY)+minY   
            absolute_maxY = ((i[3] - bottom)/(top-bottom))*(maxY-minY)+minY
            rect_locations .append((absolute_minX,absolute_minY,absolute_maxX,absolute_maxY))
    return rect_locations    

def release_rect(rect_locations,polygonName):
    k = [[0 for i in range(len(rect_locations))] for i in range(len(rect_locations))]
    
    with arcpy.da.SearchCursor(polygonName, ["SHAPE@","区域"]) as cursor:
        for row in cursor:
            maxX = max([i.X for i in row[0][0] if i ])
            maxY = max([i.Y for i in row[0][0] if i ])
            minX = min([i.X for i in row[0][0] if i ])
            minY = min([i.Y for i in row[0][0] if i ])
    
            for index,rect in enumerate(rect_locations):
                if minX > rect[0] and maxX < rect[2] and minY >rect[1] and maxY < rect[3]:
                    assert (int(row[1])-1 < len(k[index])),"图斑地域标错了"
                    k[index][int(row[1])-1] += 1    
            
    count = 0
    return_dict = {}
    for i in k:
        return_dict[ i.index(max(i))+1] = rect_locations[count]
        count += 1
    return return_dict         


#点到直线的距离，欧氏距离
def get_distance_point2line(point, line_ab):
    """
    Args:
        point: [x0, y0]
        line_ab: [k, b]
    """
    k, b = line_ab
    distance = (k * point[0] - point[1] + b) / math.sqrt(k**2 + 1)
    return distance



#程序开始点
#获取群系名称，以及地域
#str default_name = 91_祁连圆柏林_寒温带_图斑
mxd = arcpy.mapping.MapDocument('CURRENT')
lsd = mp.ListDataFrames(mxd)
default_name = ""
for df in lsd :
    for layer in mp.ListLayers(df):
        if "_图斑" in layer.name:
            assert (not default_name or default_name == layer.name),"Error:101 %s" %layer.name
            default_name = layer.name
assert default_name,"Error:102"
df = lsd[0]
#获得当前温度带，当前类型的图斑哈希
belong_tb = ""
poly_hash = {}
poly_hash2 = {}
with arcpy.da.UpdateCursor(default_name,["belong_tb","区域","类型","SHAPE@","real_TB"]) as cursor:
    for row in cursor:
        assert not(belong_tb and row[0] != belong_tb),"Error:106"
        belong_tb = row[0]
        real_tb = row[4]
        if 1:
            ll = [[i.X,i.Y] for i in row[3][0] if i]
            maxX = max([i.X for i in row[3][0] if i ])
            maxY = max([i.Y for i in row[3][0] if i ])
            minX = min([i.X for i in row[3][0] if i ])
            minY = min([i.Y for i in row[3][0] if i ])
            if  str(row[1]) in poly_hash:
                poly_hash[str(row[1])].append(ll)
                poly_hash2[str(row[1])].append([maxX,maxY,minX,minY])
            else:
                poly_hash[str(row[1])]=[ll]
                poly_hash2[str(row[1])]=[[maxX,maxY,minX,minY]]               
        
               
#dict return_dict: {1:List rect,2:List rect} Listrect:(minX,minY,maxX,maxY)
return_rect= got_rect_locations()
return_dict= release_rect(return_rect,default_name)

if target_type:
    for i in return_dict.keys()[:]:
        if i not in target_type and str(i) not in target_type:
            del return_dict[i]
        
#汇总全群系温度带图斑

savelist = []

pathList9 = []
path = poly_path.decode("utf8").encode('gbk')
def enumeration(dirname):
    for root,dirs,files in os.walk(dirname):        
        for file in files:            
            if file.endswith(".mxd") and not savelist:
                savelist.append( file.decode('gbk').encode('utf-8').split("_")[0])
                savelist.append( file.decode('gbk').encode('utf-8').split("_")[1])
            pathList9.append(os.path.join(root,file).decode('gbk').encode('utf-8'))

enumeration(path)
assert savelist,"Error 108"
savecode,savename = savelist

self_poly = []
count = 0
for i in pathList9:
    if "基础" in i and   i.endswith("带_图斑.shp") :
        count += 1
        ct2 = 0
        with arcpy.da.UpdateCursor(i,["FID",'SHAPE@',"type","name","value"]) as cursor:
            for row in cursor:
                ct2 += 1
                self_poly.append(row)
for j in ["shp","cpg","dbf","prj","sbx","sbn","shx"]:
    shutil.copyfile(total_split_temp_path.replace("shp",j),total_split_temp_path.replace("total_temp_split.shp","self_polygon.shp").replace("shp",j))
        
arcpy.mapping.AddLayer(df,arcpy.mapping.Layer(total_split_temp_path.replace("total_temp_split","self_polygon")), "TOP")       
arcpy.mapping.ListLayers(df)[0].name = "self_polygon"  #重命名    


with arcpy.da.InsertCursor("self_polygon",(["FID",'SHAPE@',"type","name","value"])) as cursor:
    for row in self_poly:
        cursor.insertRow(row)  
        
#新建图斑图层并按照范围插入图斑

now_poly = int(default_name.split("_")[0])

for i in return_dict:
    for j in ["shp","cpg","dbf","prj","sbx","sbn","shx"]:
        shutil.copyfile(total_split_temp_path.replace("shp",j),total_split_temp_path.replace("total_temp_split.shp","poly_area%s.shp"%(str(i))).replace("shp",j))
    arcpy.mapping.AddLayer(df,arcpy.mapping.Layer(total_split_temp_path.replace("total_temp_split.shp","poly_area%s.shp"%(str(i)))), "TOP")       
    arcpy.mapping.ListLayers(df)[0].name = "poly_area%s"%(str(i))  #重命名    
    #读取符合条件的图斑列表
    f = open(location_document ,'r')
    a = f.read()
    data = []
    count = 1
    for rect in a.split("\n"):
        xmax,ymax,xmin,ymin = eval(rect)
        if xmax < return_dict[i][0] or xmin > return_dict[i][2] or ymax < return_dict[i][1] or ymin > return_dict[i][3]:
            count+=1
            continue
        data.append(count) 
        count += 1
    
    poly_data = []
    with arcpy.da.UpdateCursor(total_split_path,["FID",'SHAPE@',"type","name","value","ORIG_FID"]) as cursor:
        for row in cursor:
            if int(row[-1]) in data and row[-2] != now_poly:
                for point in row[1][0]:
                    if not point:
                        continue
                    X,Y = point.X,point.Y
                    if X > return_dict[i][0] and X < return_dict[i][2] and Y > return_dict[i][1] and Y < return_dict[i][3]:
                        poly_data.append(row)
                        break
    with arcpy.da.InsertCursor("poly_area%s"%(str(i)),(["FID",'SHAPE@',"type","name","value","ORIG_FID"])) as cursor:
        for row in poly_data:
            cursor.insertRow(row)    

    split_poly = []
    with arcpy.da.SearchCursor("self_polygon",(["FID",'SHAPE@',"type","name","value"])) as cursor:
        for row in cursor:
            for point in row[1][0]:
                if not point:continue
                X,Y = point.X,point.Y
                if X > return_dict[i][0] and X < return_dict[i][2] and Y > return_dict[i][1] and Y < return_dict[i][3]:
                    split_poly.append(row)
                    break

    with arcpy.da.InsertCursor("poly_area%s"%(str(i)),(["FID",'SHAPE@',"type","name","value"])) as cursor:
        for row in split_poly:
            cursor.insertRow(row)            
            
            
            
    
                 
#新建点图层
for i in return_dict:
    for j in ["shp","cpg","dbf","prj","sbx","sbn","shx"]:
        shutil.copyfile(point_temp.replace("shp",j),point_temp.replace("temp_point.shp","area%s.shp"%(str(i))).replace("shp",j))
    
    arcpy.mapping.AddLayer(df,arcpy.mapping.Layer(point_temp.replace("temp_point.shp","area%s.shp"%(str(i)))), "TOP")       
    arcpy.mapping.ListLayers(df)[0].name = "area%s"%(str(i))  #重命名    


#按照范围插入点
field_list = ['SHAPE@', 'GRID_CODE', 'ID', 'tb', 'E_sum','p_pe', 'p1_pe1','p2_pe2', 'p3_pe3', 'p4_pe4', 'p5_pe5','p6_pe6', 'p7_pe7', 'p8_pe8', 'p9_pe9', 'p10_pe10', 'p11_pe11', 'p12_pe12', 'X', 'Y', 'p', 'pe', 'dem','t_year']
for i in return_dict:        
    with arcpy.da.InsertCursor("area%s"%(str(i)) ,field_list) as cursor:
        xmin,ymin,xmax,ymax = return_dict[i]
        xmin += 3000000
        xmax += 3000000
        xr1 = max(int(xmin/100000),4)
        xr2 = min(int(xmax/100000),50)
        yr3 = max(int(ymin/100000) ,23)
        yr4 = min(int(ymax/100000),63)

        for x in range(xr1,xr2+ 1):
            for y in range(yr3,yr4+ 1):
                df_csv = pd.read_csv(split_csv_path + "\%s,%s.csv"%(x,y))
                df_csv = df_csv[(df_csv["xx"]+3000000 >xmin)&(df_csv["xx"]+3000000 <xmax) &(df_csv["yy"]>ymin)&(df_csv["yy"]<ymax) ]
                for index, r in df_csv.iterrows():
                    row = r
                    xx = float(row[-2])
                    yx = float(row[-1])
                    l = []
                    pt = arcpy.Point(xx,yx)
                    l.append(pt)
                    for ii in row[2:]:
                        l.append(ii)
                    """if str(int(l[1])) not in ['1603','1604','1605','1606','1607','1608','1609','1610','1611','1612','1613','1614','1615','1616','1617','315','1245','1246','1420']:"""
                    cursor.insertRow(l[:-2])


              


#按范围外、置信区间外、置信区间内分为三类，0=范围外，1=置信区间内 2=范围内置信区间外
for i in return_dict: 
    xs = []
    ys = []
    if str(belong_tb) == "0":
        
        j = [csv_path + "\\" + str(belong_tb) + "." + str(i) + ".0." + str(typea) +"_" + ".".join(z for z in str(real_tb))+ ".csv" for typea in range(1,17)]
    else:
        j = [csv_path + "\\" + str(belong_tb) + "." + str(i) + ".0." + str(typea) + ".csv" for typea in range(1,17)]
    
    path = unicode(j[0], 'utf-8')
    df_temp = pd.read_csv(path)
    
    for pp in j[1:]:
        try:
            path = unicode(j, 'utf-8')
            df2 = pd.read_csv(path)
            df_temp = pd.concat([df_temp,df2])
        except:
            pass
    tb_list = list(df_temp["tb"])
    min_tb,max_tb = min(tb_list),max(tb_list)
    e_sum_list = list(df_temp["E_sum"])
    line = numpy.polyfit(tb_list,e_sum_list,1) #回归
    diss =[get_distance_point2line(point, line) for point in map(lambda x,y:(x,y) ,tb_list,e_sum_list)]
    var = numpy.var(diss)
    range_ = (-1.96*(var)**0.5,1.96*(var)**0.5) # 95置信区间
    with arcpy.da.UpdateCursor("area%s"%(str(i)),(['tb', 'E_sum',"class1"])) as cursor2:
        for row in cursor2:
            if round(min_tb,8)<=round(row[0],8)<=round(max_tb,8):
                point = row[:2]
                distance = get_distance_point2line(point, line)
                if round(min(diss),8) <= round(distance,8) < round(range_[0],8):
                    row[-1] = "2"
                elif max(min(diss),8) >= round(distance,8) > round(range_[1],8):
                    row[-1] = "2"
                elif round(range_[0],8) <= round(distance,8) <= round(range_[1],8):
                    row[-1] = "1"
                else:
                    row[-1] = "0"
                cursor2.updateRow(row)
            else:
                row[-1] = "0"
                cursor2.updateRow(row)
                
#更新比例                

for i in return_dict:
    class_dict = {}
    
    arcpy.SelectLayerByAttribute_management(default_name,"NEW_SELECTION",'"区域" = %s'%(i))
    arcpy.SelectLayerByAttribute_management("area%s"%(str(i)),"NEW_SELECTION",'"tb">= 0')
    arcpy.SelectLayerByLocation_management("area%s"%(str(i)),"COMPLETELY_WITHIN",default_name)
    with arcpy.da.UpdateCursor("area%s"%(str(i)),(['GRID_CODE','class1','SHAPE@'])) as cursor:
        for row in cursor:
            if "isin" in class_dict:  
                class_dict["isin"][int(row[1])] += 1
            else:
                class_dict["isin"]= [0,0,0,0]
                class_dict["isin"][int(row[1])] += 1
    arcpy.SelectLayerByLocation_management("area%s"%(str(i)),"COMPLETELY_WITHIN",default_name,invert_spatial_relationship = True)
    with arcpy.da.UpdateCursor("area%s"%(str(i)),(['GRID_CODE','class1','SHAPE@'])) as cursor:
        for row in cursor:
            if row[0] in class_dict:  
                class_dict[row[0]][int(row[1])] += 1
            else:
                class_dict[row[0]]= [0,0,0,0]
                class_dict[row[0]][int(row[1])] += 1
    arcpy.AddField_management("area%s"%(str(i)),'class1',"TEXT") 
    arcpy.AddField_management("area%s"%(str(i)),'ratio_0',"DOUBLE") 
    arcpy.AddField_management("area%s"%(str(i)),'ratio_1',"DOUBLE") 
    arcpy.AddField_management("area%s"%(str(i)),'ratio_2',"DOUBLE") 

    ratio_dict = {}

    for grid in class_dict:
        a0 = class_dict[grid]
        a1 = class_dict["isin"]
        a2 = class_dict["isin"]
        if str(grid) == str(now_poly):
            ratio_dict[grid] = [float(j)/float(sum(a2)) for j in [a1[0],a1[1],a1[2]]]
        else:
            ratio_dict[grid] = [float(j)/float(sum(a0+a2)) for j in [a1[0],a1[1],a1[2]]]
            
    arcpy.SelectLayerByAttribute_management("area%s"%(str(i)),"NEW_SELECTION",'"tb">= 0')
    with arcpy.da.UpdateCursor("area%s"%(str(i)),(['GRID_CODE','ratio_0','ratio_1','ratio_2',"class1","class2"])) as cursor:
        for row in cursor:
            if row[0] not in ratio_dict:
                row[5] = "1"
                cursor.updateRow(row)
                continue
            row[1:4] = ratio_dict[row[0]] 
            if row[4] == "0":
                row[5] = "0"
            elif row[4] == "1":
                if row[2] >0.6321:
                    row[5] = "1"
                elif row[2] > 0.3679:
                    row[5] = "2"
                else:
                    row[5] = "3"
            elif row[4] == "2":
                if row[3] >0.6321:
                    row[5] = "1"
                elif row[3] > 0.3679:
                    row[5] = "2"
                else:
                    row[5] = "3"
            cursor.updateRow(row)

    arcpy.SelectLayerByLocation_management("area%s"%(str(i)),"COMPLETELY_WITHIN",default_name)
    with arcpy.da.UpdateCursor("area%s"%(str(i)),(['GRID_CODE','ratio_0','ratio_1','ratio_2',"class1","class2"])) as cursor:
        for row in cursor:
            row[-1] = 1
            cursor.updateRow(row)
    
    arcpy.SelectLayerByAttribute_management("area%s"%(str(i)),"NEW_SELECTION",'"tb">= 0')
    data_d = {}
    count = 0
    with arcpy.da.UpdateCursor("area%s"%(str(i)),(['GRID_CODE','ratio_0','ratio_1','ratio_2',"class1","class2"])) as cursor:
        for row in cursor:
            count += 1
            if row[0] not in data_d:
                data_d[row[0]] = [0 for e in range(9)]
                data_d[row[0]][int(row[-2]) + 3*int(row[-1])-3] += 1
            else:
                data_d[row[0]][int(row[-2]) + 3*int(row[-1])-3] += 1
                
    cc = default_name.replace("_图斑","_地域%s_演替"%(str(i))).split("_")[2]
    dd = "带-".join(iaz for iaz in cc.split("带"))[:-1].replace("-交错","交错")
    with open(save_csv_path + "\\" + "%s_%s_%s"%(savecode,savename,dd ) +"_地域%s_演替"%(str(i))+ ".csv","a+" ) as f:
        f.write("value,c1=0&c2=1,,c1=1&c2=1,,c1=2&c2=1,,c1=0&c2=2,,c1=1&c2=2,,c1=2&c2=2,,c1=0&c2=3,,c1=1&c2=3,,c1=2&c2=3,,\n")
        for ii in data_d:
            rd = [[q,float(q)/count] for q in data_d[ii]]
            rd = sum(rd,[])
            wris = ",".join(str(q) for q in rd)
            f.write(str(ii) + "," + wris + "\n")
            
            
    arcpy.AddField_management("poly_area%s"%(str(i)),'ratio_0',"DOUBLE") 
    arcpy.AddField_management("poly_area%s"%(str(i)),'ratio_1',"DOUBLE") 
    arcpy.AddField_management("poly_area%s"%(str(i)),'ratio_2',"DOUBLE") 
    
       
    with arcpy.da.UpdateCursor("poly_area%s"%(str(i)),(['value','ratio_0','ratio_1','ratio_2',"SHAPE@","ORIG_FID"])) as cursor:
        for row in cursor:

            if str(row[0]) == str(now_poly):
                maxX = max([cc.X for cc in row[4][0] if cc ])
                maxY = max([cc.Y for cc in row[4][0] if cc ])
                if round(maxX,5) in [round(a[0],5) for a in poly_hash2[str(i)]] and round(maxY,5) in [round(a[1],5) for a in poly_hash2[str(i)]]:
                    row[1:4] = [0,1,0]
                    cursor.updateRow(row)   
                elif int(row[0]) in ratio_dict:
                    row[1:4] = ratio_dict[int(row[0])] 
                    cursor.updateRow(row)

            elif int(row[0]) not in ratio_dict:
                continue
            else:
                row[1:4] = ratio_dict[int(row[0])] 
                cursor.updateRow(row)
                